% Now start doing transition dynamics

global cfun fspace; 

npts       = 5; 

optset('chebnode', 'nodetype', 2);              % Lobatto nodes (include endpoints)


if Markuptarget < 1.1

    switch experiment
    
        case 'entry'

            fspace     = fundefn('cheb', npts, 0.75*N, 1.15*N); 

        case 'uniform'
        
            fspace     = fundefn('cheb', npts, 0.75*N, 1.15*N);    

        end
    
else
            
    switch experiment
    
        case 'entry'

            fspace     = fundefn('cheb', npts, 0.85*N, 1.15*N); 

        case 'uniform'
        
            fspace     = fundefn('cheb', npts, 0.85*N, 1.25*N);    
            
    end

end


Ngrid      = funnode(fspace); 

Mugrid     = zeros(npts, 5); 
Zgrid      = zeros(npts, 5); 
Ggrid      = zeros(npts, 5); 

nsgrid     = zeros(S, npts, 5); 

% add firms incrementally, so when we add new firms, old continue


for kt = 1 : 5

unif       = rand(Nmax, S);   % fix this once so firms incrementally enter each sector 
ns         = nssave; 

for it = 1 : npts 
    
    if Ngrid(it) > N 
        
    ms        = poissrnd(Ngrid(it) - Ngrid(it-1), S, 1);

    ns        = ns  + ms;                               % add extra firms

    z         = zsave; 

    for i = 1 : S
   
      z(ns(i) + 1 : end, i) = eps;                      % kill extra firms (set z = eps)
    
    end

    z         = sort(z, 1, 'descend');

    nsmax     = max(ns); 
    z         = z(1 : nsmax, :);                        % no reason to carry extra entries

    else
     
    ns        = nssave;     
    z         = zsave; 

    for i = 1 : S
   
      z(ns(i) + 1 : end, i) = eps;                      % kill extra firms (set z = eps)
    
    end
    
    die       = unif < (N - Ngrid(it))/N;
    
    z(die)    = eps; 
        
    z         = sort(z, 1, 'descend');
           
    ns        = sum(z > eps)';        
    
    nsmax     = max(ns); 
    z         = z(1 : nsmax, :);                        % no reason to carry extra entries
 
    end
    
    nsgrid(:, it, kt) = ns; 
    
    [Mugrid(it, kt), Zgrid(it, kt), Ggrid(it, kt)]     = solution(z, p); 
    
end
 
end

cfun           = funfitxy(fspace, kron(ones(5, 1), Ngrid), [Mugrid(:), Zgrid(:), Ggrid(:)]);

close all

nn = nodeunif(1000, Ngrid(1), Ngrid(end));

% Retrieve initial steady state before we update p.xis and p.xie

p.xis          = 0; 
p.xie          = 0; 

[~, W, L, Lp, K, B, C, Q, Z, Mu, G]                             = findequilibrium([Y; N], p);

switch experiment
    
case 'entry'

  p.xis        = 0; 
  p.xie        = -0.1105;                                          % optimal entry subsidy, leads to  0.01577% welfare gains

  optset('golden', 'tol', 1e-4); 
    
  p.xie        = golden('findxie', -0.35, 0.25, Y, N, K, C, L, options, p, cfun, fspace); 
    
  Xnew         = fsolve('findequilibrium', [Y; N], options, p);

case 'uniform' 
    
  p.xis        = Mu - 1;                                         % welfare gains of 5.58644%
  p.xie        = 0; 
   
  Xnew         = fsolve('findequilibrium', [Y*1.5; N*1.1], options, p);

end

% Find New Steady State

Ynew           = Xnew(1); 
Nnew           = Xnew(2); 

[err, Wnew, Lnew, Lpnew, Knew, Bnew, Cnew, Qnew, Znew, Munew, Gnew]  = findequilibrium(Xnew,   p);

Rnew           = p.r + p.delta; 

% Transition dynamics

xpar   = [p.beta; p.varphi; p.delta; p.psi; p.F; p.nu; p.alpha; p.theta; p.phi; p.xis; p.xie];          save xpar xpar; 

x1_ss  = [K;    N];                                                                                     save x1_ss x1_ss; 
x2_ss  = [Ynew; Cnew; Bnew; Lnew; Lpnew; Knew; Nnew; Wnew; Rnew; Qnew; Znew; Munew; Gnew];              save x2_ss x2_ss; 

dynare transition noclearall;

% Cleanup

filename = 'transition';

eval(['delete ' filename '*.m']);
eval(['delete ' filename '*.mat']);
eval(['delete ' filename '*.swp']);
eval(['delete ' filename '*.log']);
eval(['delete ' filename '*.eps']);
eval(['delete ' filename '*.pdf']);
eval(['delete ' filename '*.fig']);
eval(['delete ' filename '*.jnl']);

% Populate first period values from original steady state

Yt(1)  = Y; 
Ct(1)  = C; 
Bt(1)  = B; 
Lt(1)  = L; 
Lpt(1) = Lp; 
Kt(1)  = K; 
Nt(1)  = N; 
Wt(1)  = W; 
Rt(1)  = R; 
Zt(1)  = Z; 
Mut(1) = Mu; 

Xt     = [p.delta*K; Kt(2:end) - (1 - p.delta)*Kt(1 : end - 1)]; 

T    = length(Ct) - 1; 

Wold = 1/(1 - beta)*(log(C)  - psi/(1 + nu)*L^(1 + nu)); 

Wnew = (beta.^(0 : 1 : T - 1))*(log(Ct(2 : end)) - psi/(1 + nu)*Lt(2:end).^(1 + nu)) + beta^T/(1 - beta)*(log(Ct(end)') - psi/(1 + nu)*Lt(end)'.^(1 + nu));

dW   = (exp((1 - beta)*(Wnew - Wold)) - 1)*100;

clc

switch experiment

    case 'uniform'

        Experiment = 'Uniform Subsidy';

    case 'entry'

        Experiment = 'Entry Subsidy';

end

fprintf('<strong> Table 7: Implications of Alternative Policies, Oligopoly </strong> \n'); 

fprintf('\n');

fprintf(' Steady state comparisons in percent - aggregate markup at %1.2f \n', Mu);  

fprintf('\n');

fprintf([' Experiment --> ' Experiment, ' \n']);   

fprintf('\n');

disp(table(round((Ynew/Y   - 1)*100, 1), round((Cnew/C   - 1)*100, 1), round((Lnew/L   - 1)*100, 1), round((Nnew/N   - 1)*100, 1), round((Knew/K   - 1)*100, 1), round((Znew/Z   - 1)*100, 1), round(dW, 2), ...
    'VariableNames', {'Y','C','L','N','K','Z','Welfare,%'}))

